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ABSTRACT 



Context. The migration of growing protoplanets depends on the thermodynamics of the ambient disc. Standard modelling, using 
locally isothermal discs, indicate in the low planet mass regime an inward (type-I) migration. Taking into account non-isothermal 
effects, recent studies have shown that the direction of the type-I migration can change from inward to outward. 
Aims. In this paper we extend previous two-dimensional studies, and investigate the planet-disc interaction in viscous, radiative discs 
using fully three-dimensional radiation hydrodynamical simulations of protoplanetary accretion discs with embedded planets, for a 
range of planetary masses. 

Methods. We use an explicit three-dimensional (3D) hydrodynamical code NIRVANA that includes full tensor viscosity. We have added 
implicit radiation transport in the flux-limited diffusion approximation, and to speed up the simulations significantly we have newly 
adapted and implemented the FARGO-algorithm in a 3D context. 

Results. First, we present results of test simulations that demonstrate the accuracy of the newly implemented FARGO-method in 3D. 
For a planet mass of 20 Memth we then show that the inclusion of radiative effects yields a torque reversal also in full 3D. For the same 
opacity law used the effect is even stronger in 3D than in the corresponding 2D simulations, due to a slightly thinner disc. Finally, we 
demonstrate the extent of the torque reversal by calculating a sequence of planet masses. 

Conclusions. Through full 3D simulations of embedded planets in viscous, radiative discs we confirm that the migration can be 
directed outwards up to planet masses of about 33 Memth- Hence, the effect may help to resolve the problem of too rapid inward 
migration of planets during their type-I phase. 

Key words, accretion discs - planet formation ~ hydrodynamics - radiative transport 



1. Introduction 

The process of migration in protoplanetary discs allows form- 
ing planets to move away from their location of creation and 
finally end up at a different position. The cause of this change 
in distance from the star are the tidal torques acting from the 
disturbed disc back on the protoplanet. These can be separated 
into two parts: i) the so-called Lindblad torques that are cre- 
ated by the two spiral arms in the disc and ii) the corotation 
torques that are caused by the co-orbital material as it periodi- 
cally exchanges angular momentum with the planet on its horse- 
shoe orb it. For comprehensive in tro ductions t o the field see for 
example iPapaloizou et alj ( l2007h or iMassetl ( l2008h . and refer- 
ences therein. The Lindblad torques, caused by density waves 
launched at Lindblad resonances, quite generally lead to an in- 
ward motion of the planet explaining q uite nicely the existence 
of the observed hot planets dWardll 19971) . The corotation torques 
are caused mainly by two e ffects, first by a gradient in the 
vortensity (Tanaka et al. 2002), and second by a gradient in the 
entropy CBaruteau & Masset 2008i) . For typical protoplanetary 
discs both contributions can be posi tive, possibly counterbal- 
ancing the negative Lindblad torques dPaardekooper & Mellemal 
[2OO6; Baruteau & Masset 2008). For the typically considered lo- 
cally isothermal discs where the temperature depends only on 
the radial distance from the star the net torque is ne gative and mi- 
gratio n directed inwards for typical disc parameter dTanaka et alj 
I2OO2I) . Through planetary synthesis models the inferred rapid in- 



ward migration of planetary cores has been found to be inconsis- 
tent with the observed mass-distanc e distribution of exoplanets 
dAlibert et al.ll2004l: llda"& Linll2008[) . Possible remedies are the 
retention of icy cores at the snow line or a strong reduction in 
the speed of type-1 migration (of embedded small mass planets). 
Here we focus of the latter process. 

Different mechanism for slowi ng down the too rapid in- 
ward migration have been discussed (jMas set et al.l2006tlLi et al.l 
,2009) , but the inclusion of more realistic physics seems to 
be particularly appeaUng. As h as been pointed out first by 
iPaardekooper & Mellemal d2006l) the inclusion of radiative trans- 
fer can cause a strong reduction in the migration speed. 
The pro cess has been subsequent l y investigated by several 



group s ( Baruteau & Massei 2008 
I2OO8I : IPaardekooper & Papaloizottl 



Paardekooper & Mellemal 



2008t iKlev & Cridal l2008h . 
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who show that indeed the migration process can be slowed 
down or even reversed for sufficiently low mass planets. This 
new effect occurs in non-i sothermal discs and scal es with 
the gradient of the entropy dBaruteau & Massed l2008l) . hence 
entropy-related torque. However, in a strictly adiabatic situ- 
ation after a few libration time scales the entropy gradient 
will flatten within the corotation region due to phase mix- 
ing. This will lead to a saturation (and subsequent disappear- 
ance) of the part of the corotation torque t hat is caused by the 
entropy gradient in the horseshoe regio n dBaruteau & Masseti 
I2008t IPaardekooper & Papaloizoull200"9l) . To prevent saturation 
of this entropy-related torque som e radiative diffusion (or lo- 
cal radiative cooling) is required dKlev & Cridal I2OO8I) . Since 
for genuinely inviscid flows, the streamlines in the horse- 
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shoe region will be closed and symmetric with respect to the 
planet's location, some level of viscosity is always necessary to 
avoid torque saturation ( Masset 2001; Ogilvie & Lubow 20031 
iPaardekooper & Papaloizoull2008i 120091) . This applies to both, 
the vortensity- and entropy-related corotation torques. The max- 
imum planet mass for which a change of migration may oc- 
cur due to this effect lies for typical disc masses in the range 
of about 40 earth masses, beyond which the migration rate 
follows the standard (isothermal) case, as gap formation set s 
in which reduces the corotation effects jKlev & Cridal |2008|) . 
Most of the above simulations have studied only the two- 
dimensional case, while three dimensional models including ra- 
diative effects have been present ed only for very small masses 
( Paardekooper & Mellemal l2006l) . or for Jupiter type planets 
(iKlahr & Klevll2006h . A range of planet masses has not yet been 
studied systematically in full 3D. 

In this paper we investigate the planet-disc interaction in ra- 
diative discs using fully three-dimensional radiation hydrody- 
namical simulations of protoplanetary accretion discs with em- 
bedded planets for a variety of planetary masses. For that pur- 
pose we modified and substantially extended a n existing multi- 
dimen s ional hydrodyna mical code Nirvana dZiegler & Yorkd 
Il997al; iKlev et al.ll200lb by incorporating the FARGO-algorithm 
dMassetl 12000^ " and ra diative transport in the flux-l i mited 
diffus ion approximation (iLevermore & Pomraningiri98ll: iKlevI 
119891) . The code Nirvana can in principle handle nested grids 
which allows to z oom-in on the detailed struc ture in the vicin- 
ity of the planet jP'Angelo et al.ll2002l l2003b . however in the 
present context we limit ourselves to single grid simulations. We 
present several test cases to demonstrate first the accuracy of the 
FARGO-method in 3D. We then proceed to analyse the effects of 
radiative transport on the disc structure and torque balance. For 
our standard planet of 20Meanh, we find that the effect of torque 
reversal appears to be even stronger in 3D than in 2D for an oth- 
erwise identical physical setup. We have a more detailed look at 
the implementation of the planet potential and show that it has 
definitely an influence on the strength of the effect. Finally, we 
perform simulations for a sequence of different planet masses 
to evaluate the mass range over which the migration may be re- 
versed. The consequence for the migration process and the over- 
all evolution of planets in discs is discussed. 



2.1. Basic equations 

Discs with embedded planets have mostly been modelled 
through 2D simulations in which the disc is assumed to be in- 
finitesimal thin, and vertical integrated quantities are used to de- 
scribe the time evolution of the disc with the embedded planet. 
This procedure saves considerable computational effort but is 
naturally not as accurate as truly 3D simulations, in particular 
the radiation transport is difficult to model in a 2D context. 

In this work we present an efficient method for 3D disc sim- 
ulations based on the FARGO algorithm (lMassetll2000"i) . For ac- 
cretion discs where material is orbiting a central object the best 
suited coordinates are spherical polar coordinates (r, 6, if) where 
r denotes the radial distance from the origin, 9 the polar angle 
measured from the z-axis, and ip denotes the azimuthal coordi- 
nate starting from the jc-axis. 

In this coordinate system, the mid-plane of the disc coincides 
with the equator {9 - njl), and the origin of the coordinate sys- 
tem is centred on the star. Sometimes we will need the radial 
distance from the polar axis which we denote by a lower case s, 
which is the radial coordinate in cylindrical coordinates. 

For a better resolution of the flow in the vicinity of the planet, 
we work in a rotating coordinate system which rotates with the 
orbital angular velocity Q, which is identical to the orbital angu- 
lar velocity of the planet 



G (M, + nip) 



1/2 



(1) 



where M« is the mass of the star, nip the mass of the planet, and 
a the semi-major axis of the planet. Only for testing purposes for 
our implementation of the FARGO-method we let the planet move 
under the action of the disc. 

The Navier-Stokes equations in a rotating coordinate system 
in spherical coordinates read explicitly: 
a) Continuity Equation 



dp 



(2) 



Here p denotes the density of the gas and u = (uy, ug, u^) its 
velocity. 

b) Radial Momentum 



2. Physical modelling 

The protoplanetary disc is treated as a three-dimensional 
(3D), non-self-gravitating gas whose motion is described 
by the Navier-Stokes equations. The turbulence in discs is 
thou ght to be driven by m agneto-hydrodynamical instabili- 
ties (iBalbus & Hawlevlll998h . Since we are interested in this 
study primarily on the average effect the disc has on the 
planet, we prefer in this work to simplify and treat the disc 
as a viscous medium. The dissipative effects can then be de- 
scribed via the standard viscou s stress-tensor approach (eg. 
iMihalas & Weibel MihalasTT98l . We assume that the heating of 
the disc occurs solely through internal viscous dissipation and 
ignore in the present study the influence of additional energy 
sources such as irradiation from the central star or other external 
sources. The internally produced energy is then radiatively dif- 
fused through the disc and eventually emitted from its surfaces. 
To describe this pro cess we utihse th e flux-limited diffusion 
approximation (FLD, ILevermore & Po mraning 198 li) which al- 
lows to treat approximately the transition from optically thick to 
thin regions near the disc's surface. 



dpur 



+ V-{pUr\x) - p — + pr sin 6* (w + Q) 
ot r 

dp ^ 
+ pa, - — -p—+f,. 
or or 



(3) 



Here co is the azimuthal angular velocity as measured in the 
rotating frame, p is the gas pressure, and O denotes the grav- 
itational potential due to the star and the planet. The vector 
a - [a,; ag, a^) represents inertial forces generated by the ac- 
celerated origin of the coordinate system. Specifically, a equals 
the negative acceleration acting on the star due to the planet(s). 
c) Meridional Momentum 



d prug 
dt 



+ V ■ iprugvi) - pr" sin COS 0(w + Q) 



dp 50) 

+ prog — p — — ^ r fg 

^ 89 ^ 89 



(4) 



d) Angular Momentum 
8ph, 



dp 8^ 



+ V ■(ph,u) =prsin9a^- p— +rsm9L; (5) 

dip dip 
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where we defined the total specific angular momentum (in the 
inertial frame) 



h, = r sin (w + Q) ; 



(6) 



i.e. the azimuthal velocity in the rotating frame is given by - 
(jL>r sin0. 

The Coriolis force in Eq. (|5]) for (or h,) has been incorpo- 
rated into the left hand side. Thus, it is written in such a way as to 
conserve total angular momentum best. This conservative treat- 
ment is necessar y to obtain an accurate solution of the embedded 
planet problem (lKlevlll998l) . 

The function f = (fr,fe,fip) in the momentum equations 
denotes the viscous forces which are stated explicitly for the 
three-d imensional case in spherical polar coordinates in iTassoull 
(Il978l) . For the description of the viscosity we use a constant 
kinematic viscosity coefficient v. 
e) Energy equation (internal energy) 



dpCyT 
dt 



+ V ■ {pcyTn) = -pV ■ 



(7) 



Here T denotes the gas temperature in the disc and Cy is the spe- 
cific heat at constant volume. On the right hand side, the first 
term describes compressional heating, the viscous dissipa- 
tion function, and F denotes the radiative flux. In writing eq. (|7]i 
we have assumed that the radiation energy density E - urT'^ 
is always small in comparison to the thermal energy density 
e - pc^J . Here, denotes the radiation constant. Furthermore, 
we utilise the one-temperature approach and write for the radia- 
tive flux, using flux-limited diffusion (FLD) 



F ^ —VT, 

P(k + o") 



(8) 



where c is the speed of light, k the Rosseland mean opac- 
ity, cr the scattering coefficient, and A the flux-limiter Using 
FLD allows us to perform stable accretion disc models that 
cover several vertica l scale heights. We use h e re the FLD ap- 
proach descri bed in 
flux-limiter of iKlevI d 



Levermore & Pomranind (Il98lb with the 



1989|). Its suitabilitv for protostellar discs 



has been shown in iKlev&LinI (1 996. 199^, and for embed- 
ded planets in iKlahr & KlevH20b6l) . In this work we use for 
the Rosseland mean opa city k(p, T) the analytical formulae by 
lLin&Papaloizoul(ll985h and set the scattering coefficient cr to 
zero. To close the basic system of equations we use an ideal 
gas equation of state where the gas pressure is given hy p - 
RgaspT/fi, with the mean molecular weight yu and gas constant 
Rgas- For a standard solar mixture we assume here yu = 2.35. The 
speed of sound is calculated from Cs - yfypTp with the adiabatic 
index y - 1.43. 



2.2. Planetary potential 

The total potential 3) acting on the disc consists of two contribu- 
tions, one from the star O,, the other from the planet Op 



Gnir, 



p 



■rp)2 



where Tp denotes the radius vector of the planet location. The 
embedded planet is modelled as a point mass that orbits the cen- 
tral star on a fixed circular orbit. In numerical simulations the 
planetary potential has to be smoothed over a few gridcells to 
avoid divergences. 
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Fig. 1. The gravitational potential of a 20Meaiih planet with dif- 
ferent smoothings applied, see Eqs. (|9]l and ( fTOl i. The distance 
d from the planet is given in units of Op (here ajup), and the 
smoothing length r,,,, in units of the Hill radius of the planet 
which refers here to Rh = 0.027 la^. Note that the data points 
are drawn directly from the used 3D computational grid and in- 
dicate our standard numerical resolution (see section 12741 . 



Typically in 2D simulations the planetary potential is mod- 
elled by an e-potential 



Of, = — 



Gm„ 



(9) 



where we denote the distance of the disc element to the planet 
with d = |r - r,,| and e is the smoothing length. In a 2D con- 
figuration a potential of this form is indeed very convenient, as 
the smoothing takes effects of the otherwise neglected vertical 
extent of the disc into account. To account for the finite disc 
thickness H which depends on the temperature in the disc, an 
often used value for the smoothing length in 2D-simulations is 
e = fsm = 0.6//. The obtained 2D torques are then found to be 
in reasonable agreement with three-dimensional analytical esti- 
mates of the torques, that do not include a softening length for 
the planet potential. 

However in a 3D configuration the same approach is not 
necessary and would lead to an unphysical 'spreading' of the 
potential over a large region. Hence, we ap ply a different type 
of smoothing, follow iKlahr & KlevI (|2006|) and use a cubic- 
potential of the form 



O 



,cub 



d 



d 



for d < Ts, 
for d > r.. 



(10) 



This potential is constructed in such a way as to yield for dis- 
tances d larger than rs,n the correct 1 /r potential of the planet, 
and inside that radius {d < r^m) the potential is smoothed with 
a cubic polynomial such that at the transition radius r^m the po- 
tential and its first and second derivative agree with the analytic 
outside 1/r-potential. To illustrate the various types of smooth- 
ing, we display in Fig. [1] the behaviour of the different forms 
of the planetary potential. Clearly the e-potential leads for the 
same values of r^m to a much wider and shallower potential than 
our cM/j/c-approach. For the often used value e - 0.6// the e- 
potential is felt way outside the Hill radius 



R^ 



3M. 



1/3 
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and leads to a significant underestimate of the potential depth al- 
ready at Tsni- The cM/j/c-potential (Eq.fTOji will always be accurate 
down to li = and is inside much deeper than the e-potential, 
and hence more accurate. In the simulations presented below we 
study in detail the influence that the potential description has on 
the value of the torques acting on the planet. 

We calculate the gravitational torques acting on the planet by 
integrating over the whole disc, where we apply a tapering func- 
tion to exclude the inner parts of the Hill sphere of the planet. 
Specifically, we use the smooth (Fermi-type) function 



Md) = 



exp 



d/Rn-b 
b/lO 



+ 1 



(11) 



which increases from at the planet location (d = 0) to 1 outside 
d > Rh with a midpoint fb - I /2 at d = bRn, i.e. the quantity b 
denotes the torque-cutoff radius in units of the Hill radius. This 
torque cutoff is necessary to avoid first a large, possibly very 
noisy contribution from the inner parts of the Roche lobe, and 
second to disregard material that is gravitationally bound to the 
planet. The question of torque cutoff and exclusion of Roche 
lobe material becomes very important when (/) the disc self- 
gravity is neglected, and (;;) there exists material bound to the 
planet (e.g. a circumplanetary disc). This issue should definitely 
be addressed in the futu re. Here, we assu me a transition radius 
of ^ = 0.8 Hifl radii (see lCrida et aLll2008l Fig. 2). For reference 
we quote the width of the horseshoe region which is given for an 
isothermal disc approximately by dMasset et al.ii200q) 



1.16fl 



(H/r) 



(12) 



with the mass ratio q = nip/Mt and the local relative disc thick- 
ness H/r. For an adiabat ic disc H has to be replaced by jH 
dBaruteau & Massetll2008l) . 



2.3. Setup 

The three-dimensional (r, 6, (f) computational domain consists of 
a complete annulus of the protoplanetary disc centred on the star, 
extending from rmm - 0.4 to rmax = 2.5 in units of ro - ajup - 
5.2AU. In the vertical direction the annulus extends from the 
disc's midplane (at 9 = 90°) to about 1° (or 6 = 83°) above 
the midplane. In case of an inclined planet the domain has to 
be extended and cover the upper and lower half of the disc. The 
mass of the central star is one solar mass M, = M©, and the total 
disc mass inside [r,nin, r^ax] is M^isc - O.OIM©. For the present 
study, we use a constant kinematic viscosity coefficient with a 
value of V = lO'^cm^/s, a value that relates to an equivalent 
a - 0.004 at ro for a disc aspect ratio of H/r = 0.05, where 
V - aH^Q.K. In standard dimensionless units we have v — 10"^. 

The models are initialised with a locally isothermal config- 
uration where the temperature is constant on cylinders and has 
the profile r(i) oc i ' , where s is related to r through i = r sin 9. 
This yields a constant ratio of the disc's vertical height H to 
the radius s. The initial vertical density stratification is approxi- 
mately given by a Gaussian: 



p{r, 9) = Pair) exp 



{n/2 - 9Y r^ 



//2 



(13) 



Here, the density in the midplane is po('") °^ ''^ which leads 
to a S(r) oc r^^l^ profile of the vertically integrated surface 
density. The vertical and radial velocities, ug and m,-, are ini- 
tialised to zero. The initial azimuthal velocity is given by 
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Fig. 2. Evolution of semi-major axis, eccentricity and inclination 
as a function of time for a 20 Mg^rih planet in a three-dimensional 
disc. Results are displayed for different numerical resolutions 
(Nr, Ng, N^), and using Fargo and Non-Fargo runs. 



the equilibrium of gravity, centrifugal acceleration and the ra- 
dial pressure gradient. In case of purely isothermal calculations 
this setup is equal to the equilibrium configuration (in the case 
of closed radial boundaries). For fully radiative simulations the 
model is first run in a 2D axisymmetric mode to obtain a new 
self-consistent equilibrium where viscous heating balances ra- 
diative transport/cooling from the surfaces (see Sect. |4.1| below). 
This initialisation through an axisymmetric 2D phase (in the r-9 
plane) reduces the required computational effort substantially, as 
the evolution from the initial isothermal state towards the radia- 
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live equilibrium takes about 100 orbits, if the disc is started with 
an isothermal equilibrium having constant Hjr. After reaching 
the equilibrium between viscous heating and radiative trans- 
port/cooling, we extend this model to a full 3D simulation, by 
expanding the grid into the ^-direction, and the planet is embed- 
ded. 

2.4. Numerics 

We adopt a coordinate system, which rotates at the orbital 
frequency of the planet. For our standard cases, we use an 
equidistant grid in r,6,(p with a resolution of {Nr,Ne,N^) - 
(266, 32, 768). To minimise disturbances (wave reflections) from 
the radial boundaries, we impose, at r^in and rmax, damping 
boundary conditions where all three velocity components are re- 
laxed towards their initial state on a timescale of approximately 
the local orbital period. The radial velocities at the inner and 
outer radius vanish. The angular velocity is relaxed towards the 
Keplerian values. For the density and temperature, we apply re- 
flective radial boundary conditions. In the azimuthal direction, 
periodic boundary conditions are imposed for all variables. In 
the vertical direction we apply outflow boundary conditions. The 
boundary conditions do not allow for mass accretion through the 
disc, such that the total disc mass remains nearly constant dur- 
ing the time evolution, despite a possible small change due to 
little outflow through the vertical boundaries and the used den- 
sity floor (see below). 

The numerical details of the used finite volume code 
(NIRVANA ) relevant for thes e pl anet disc simulations w ere de- 
scribed in lKlev et all (1200 ll) and lD'Angelo et al.l ( |2003l) . In the 
latter paper the usage of the nested grid-technique is described in 
more detail as well. The original version of the NIRVANA code, 
on which our programme is based upon, has been developed 
bv IZiegler & Yorke (1997b). The empowerm ent with FARGO is 
based on the original work bv lMassel (l2000al) . Our implementa- 
tion appears to be the first inclusion of the FARGO-algorithm in 
a 3D spherical coordinate system. More details about the imple- 
mentation are given in the appendix. The basic algorithm of the 
newly implemented radiation part in the energy equation (|7]i is 
presented in the appendix as well. To avoid possible time step 
limitations this part is always solved implicitly. 

3. Test calculations 

3.1. The FARGO -a /gor/fhm 

To test the 3D implementation of the FARGO-algorithm in our 
NIRVANA-code we have run several models with planets on cir- 
cular, elliptic and inclined orbits with and without the FARGO- 
m ethod applied. Here, we follow closely the models presented 
in ICresswell et al.l ( |2007|) and consider moving planets in 3D 
discs. As the tests are dynamically already complicated we use 
here only the isothermal setup. The different setups gave very 
similar results in all cases, and we present results for one com- 
bined case of a 20 Mgaith planet embedded in a locally isother- 
mal disc with an initial non-zero eccentricity (e - 0.2) and 
non-zero inclination (/ = 5°). All phy sical parameters of thi s 
run are identical to those described in ICresswell et al.l (|2007|) . 
and we compare our results to the last models presented in 
that paper (their Fig. 16). The outcome of this comparison is 
shown in Fig.|2l where we display the results of a standard non- 
Fargo run with the r esolution (Nr, No , N^) - (264, 80, 800) with 
the data taken from ICresswell et alj (|2007|) (where a different 
code has been used) to two runs having a lower resolution of 
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Fig. 3. Radial stratification of surface density (left) and the mid- 
plane temperature (right) in the disc. This dashed lines repre- 
sent simple approximations to the 3D stratified results. The solid 
(green) curves labelled 'flat' refer to corresponding results for a 
vertically integrated flat 2D disc using the same input physics. 
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Fig. 4. Vertical stratification of density (left) and temperature 
(right) at a radius of r = 1.44. Thin dashed lines just represent 
simple functional relations. 



{Nr, No, N^) = (128, 34, 384), one with FARGO and the other one 
without. We can see that all 3 models (obtained with two differ- 
ent codes, methods and resolutions) yield very similar results. 
The scatter of the data points is slightly smaller in the FARGO- 
run. 

3.2. The radiation algoritiim 

To obtain an independent test of the newly implemented radia- 
tion transport module in our NIRVANA-code we performed a run 
with the standard setup as described above but with no embed- 
ded planet. Hence, this setup refers to an axisymmetric disc with 
internal heating and radiative cooling. For a fixed, closed compu- 
tational domain it is only the total mass enclosed that determines 
the final equilibrium state of the system, once the physics (vis- 
cosity, opacity, and equation of state) have been prescribed. The 
radial dependence of the vertically integrated surface density and 
the midplane temperature are displayed in Fig.|3] and the corre- 
sponding vertical profiles at a radius of r = 1.44 in Fig. |4| First 
of all, these new results obtained with NIRVANA agree very well 
with those obtained with the completely independent 2D-code 
RH2D used in the r - 6 mode, such as presented for example in 
iKlev et al] (Il993h . which are not shown in the figures, however 
As the final configuration of the system is given by the equilib- 
rium of internal (viscous) dissipation and radiative transport, this 
test demonstrates the consistency of our implementations. 

To relate our 3D results to previous radiative 2D runs which 
use vertically integrated quantities, and hence can only use 
an ap proximative energy transport and cooling dKlev & Cridal 
l2008l) . we compare in Fig. [3] the results obtained with the two 
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methods. Both models are constructed for the same disc mass 
and identical physics. The label 'flat' in the Figure refers to the 
flat 2D case (obtained with RH2D, see Kley & Crida (2008|)) and 
the 'stratified' label to our new 3D implementation presented 
here. The left graph displays the vertically integrated surface 

density distribution, here = pr sm(6)d6, with /i = 1 

for two-sided and 2 for one-sided discs. Our result is well repre- 
sented by a S oc r '^^ profile as expected for a closed domain and 
constant viscosity. Interesting is the irregular structure at radii 
smaller than r x 1 .0 in the full 3D stratified case, and we point 
out that these refer to the onset of convection inside that radius. 
To model convection is of course not possible in a flat 2D ap- 
proach. The temperature distribution for the full 3D case follows 
approximately a T oc r profile. Here, the approximate flat- 
disc model leads to midplane temperatures that are about 40% 
higher for the bulk part of the domain than in the true 3D case. 
Possibly a refined modelling of the vertical averaging procedure 
and the radiative losses in the flat 2D case could improve the 
agreement here, but in the presence of convection we may ex- 
pect differences in any case. 

In Fig.|4]we display the vertical stratification of the disc at 
a specific radius in the middle of the computational domain at 
r - 1 .44. Two simple approximations are over-plotted as dashed 
lines. Note, that in these plots the stratification is plotted along 
the r = const, lines which deviates for thin discs only slightly 
from z - rcosO. Taking zq - 0.08, the Gaussian curve for 
the density refers to po exp [-(z/zo)"] and the temperature fit to 

T(z) = ro[l - 0.4(z/zo)^]- These simple formulae are intended 
to guide the eye rather than meant to model exactly the struc- 
ture at this radius which depends on the used opacities. Given 
the simplicity of these, it is interesting that they approximate the 
true solution reasonably well within one scale height. 

3.3. Density Floor 

By expanding the computed area in the 0-direction beyond the 
90° to 83° region of our standard model, the code would have 
to cover several orders of magnitude in the density. Thus, many 
more grid cells would be required to resolve the physical quanti- 
ties. In order to avoid this and save computation time, we apply 
a minimum density function (floor) for the low-density regions 
high above the equatorial plane of the disc. It reads 



P 

Pmin 



for p > p„„„ 
for p < p„„„ 



(14) 



Of course, applying a density floor like this will create mass 
inside the computed domain. The density floor p,„„, has now to 
be chosen such that: firstly the computation is not handicapped 
by too small values and secondly the inner (optically thick) parts 
of the disc are not influenced. To test the sensitivity of the disc 
structure against the density floor we performed a series of test 
calculations, and show the results of simulations with diff'erent 
p,„i„ in Fig.|5] These runs cover - 90° to 70°, a range about 3 
times as large as before. The density and temperature profiles in 
these simulations do not differ for the regions near the equatorial 
plane, because the density is too high for the minimum density 
to take effect. Indeed, all curves are nearly indistinguishable in 
the region for optical depths larger than t - 1.0, with 



pKdz. 



(15) 
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Please note, that in the plot we do display the results along lines 
of constant (spherical) radius. Moving further away from the 
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Fig. 5. Vertical stratification of density in logarithmic scale (top) 
and temperature (bottom) at a radius of r = 1 .00 for a simulation 
covering the - 90° to 70° region. The optical depth t - 1 .0 is 
reached at about z - 0.055. 



equatorial plane, one can see in the density profile the different 
minimum densities, but in the temperature profile there is hardly 
any difference at all. Above a certain distance from the equatorial 
plane the temperature remains constant. The little fluctuations 
visible in the profile are due to oscillations in the temperature 
for the low mass regions. 

By applying a minimum density the code is capable of re- 
solving large distances above the equatorial plane with a reason- 
able number of grid cells. Also note that it is not necessary to use 
a minimum density for calculations covering only the 6 - 90° to 
83° regions, as the density is always high enough. 



4. Models with an embedded planet 

For all the models with embedded planets we use our standard 
disc setup as described in section 12.31 with the corresponding 
boundary conditions in section [Z4l Here, we briefly summarise 
some important parameter of the setup. The three-dimensional 
(r, 0, (fi) structure of the disc extends form rmin - 0.4 to r,nax =2.5 
in units of ro = ay,,^ = 5.2AU. In the vertical direction the an- 
nulus extends from the disc's midplane (at = 90°) to about 7° 
(or - 83°) above the midplane. For our chosen grid size of 
{Nr,Ncp,No) = 266,32,768) this refers to linear grid resolution 
of A a! 0.008 at the location of the planet, which corresponds to 
3.3 gridcells per Hill radius, and to about 5 gridcells per horse- 
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Fig. 6. Density (top) and Temperature (bottom) for a fully radia- 
tive model in a 2D axisymmetric simulation 



shoe half-width for a 20Meanh planet in a disc with H/r - 0.05. 
In this configuration the planet is located exactly at the corner of 
a gridcell. In the fully radiative disc, the temperature at the disc 
surface is kept at the fixed ambient temperature of 10 K. This 
simple Tow-temperature' boundary condition ensures that all the 
internally generated energy is liberated freely at the disc's sur- 
face. It is only suitable for optically thin boundaries and does not 
influence the inner parts of the optically thick disc (see Fig. |5]l. 
The disc has a mass of 0.0 IM©, and an aspect ratio Hjr - 0.05 
in the beginning. 

4.1. Initial setup 

Before placing the planet into the 3D disc we have to bring it 
first into a radiative equilibrium state such that our results are 
not corrupted by initial transients. As described above this ini- 
tial equilibration is performed in an axisymmetric 2D setup that 
is then expanded to full 3D. Tests with our code have shown that 
we reach the 3D equilibrium state (a constant torque) in a cal- 
culation with embedded planets about 50% faster when starting 
first with the 2D radiative equilibrium disc. 

In Fig. |6] the 2D density and temperature distributions for 
such an equilibrated disc are displayed. In the equilibrium state 
of the fully radiative model the disc is much thinner in compar- 
ison to the isothermal starting case, see Fig. |2l Consequently, 
the density is increased in the equatorial plane, leaving the areas 
high above and below the disc with less material. Apparently, 
for this disc mass and the chosen values of viscosity and opacity, 
the balance of viscous heating and radiative cooling reduces the 
aspect ratio of the disc from initially 0.05 to about 0.037 in the 
radiative case. Had we started with an initially thinner disc, the 
difference would of course not be that pronounced. 

After successfully completing the equilibration we now em- 
bed a 20 Mearth planet into the disc. The planet is held on a fixed 
orbit and we calculate the torques acting on it through integrating 
over the whole disc taking into account the above tapering func- 
tion with a cutoff riorq - 0.8/?//, which refers tob - 0.8 in Eq.fTTl 
In addition to this value of the torque cutoff he have tested how 
the obtained total torque changes when using b - 0.6. For our 
standard 20Meaiih planet presented in the following we found that 
for the isothermal cases the results change by about 10% and in 
the radiative case by about 30%, which can be considered as 
a rough estimate of the numerical uncertainties of the results. 
The deviation is larger in the radiative situation because in this 
case important (corotation) contributions to the total torque orig- 
inate from a region very close to the planet, which is influenced 
stronger by the applied torque cutoff. Here, cancellation effects 
caused by adding the negative Lindblad and the positive corota- 
tion torque may explain part of the larger relative uncertainty in 
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Fig. 7. Vertical density distribution at r = 1 .4 for a fully radiative 
model in a 2D axisymmetric simulation, 'h-': isothermal starting 
configuration, : relaxed radiation equilibrium configuration 



the radiative case. We note, that our applied torque cutoff" is not 
hard but refers to the smooth function (fTTT i. Keeping in mind that 
there are only 3.3 gridcells per Hill radius, smaller values for b 
are not useful. 



4.2. Isothermal discs 

Due to the applied smoothing, we expect the planetary potential 
to modify the density structure of the disc near the planet and 
subsequently change the torques acting on the planet. First, we 
investigate the influence of the planetary potentials (see Fig. [1]) 
on the disc and torques in the isothermal regime. The 2D surface 
density distribution in the disc's midplane at 100 planetary or- 
bits corresponding to our two extreme planetary potentials (the 
shallowest and the deepest) is displayed in Fig. [8] where we used 
a cutoff for the maximum displayed density to make both cases 
comparable. As expected, a deeper planetary potential results in 
a higher density concentration inside the planetary Roche lobe 
and to a slightly reduced density in the immediate surroundings. 
This accumulation of mass near the planet for deeper potentials 
is illustrated in more detail Fig. |9] For our deepest r,,,, = 0.5 
cubic potential the maximum density inside the planet's Roche- 
lobe is over an order of magnitude larger than in the shallowest 
rsm = 0.8 e-potential. 

In Fig.[TO]we show the specific torques acting onto the planet 
using different potentials for the case of H/r - 0.05. The total 
torque is continuously monitored and plotted versus time in the 
upper panel. The radial torque density r(r) for the same models 
is displayed in the lower panel. Here, r(r) is defined such that 
the total torque T"" acting on the planet is given by 



r(r)dr. 



(16) 



The time evolution of the total torque displays a characteristic 
behaviour Starting from the axisymmetric case, a first interme- 
diate plateau is reached at early times between r a; 5 - 10, after 
which the torques settle on longer timescales towards their final 
equilibria. The initial plateaus correspond to the values of the 
torques shortly after the disc material has started its horseshoe- 
type motion in the co-orbital region. The level of this so-called 
unsaturated torque depends on the local disc properties and on 
the applied smoothing of the potential, as indicated clearly in the 
top panel of Fig. [10] In the following evolution, the material in 
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Fig. 8. Surface density distribution for isothermal simulations 
with HIr - 0.05 at 100 planetary orbits. Displayed are results 
for the shallowest and deepest potential. Top; e-potential with 
Vsm - 0.8, Bottom: cubic with rj,,, - 0.5. 



where Pp is the orbital period of the planet, and Xs the half- 
width of the horseshoe region (see eq. [T2l i. In our case (for 
^ = 6 X 10"^ and HIr - 0.05) the libration time is about 30P. 
The different initial values of the unsaturated torques depend 
on the form of the potential (eg. smoothing length), but note 
that the timescale to reach equilibrium is similar in all cases. 
This particular time behaviour of the torques and the process 
of saturation has been describe d recen tly for isot hermal discs 
by Paardekooper & Papaloizoul ( l2009h . see also iMasset et~aLl 
(l2006l) ^ 

The two runs with the e-potential result in the most nega- 
tive torque values, i.e. the fastest inward migration (lower two 
curves in the upper panel). While the total torques of the two e- 
potentials are nearly identical, in the corresponding radial torque 
distribution the cases are clearly separated, a fact which is due 
to cancellation effects when adding the inner (positive) and outer 
(negative) contribution. The slightly deeper cubic r,,,, = 0.8 po- 
tential leads to a marginally decreased (in magnitude) torque 
compared to the simulations with e-potential. For the cubic 
Vsm - 0.5 potential we obtain an even less negative equilibrium 
torque compared to all the other isothermal simulations. As most 
of the corotation torque is generated in the vicinity of the planet, 
a change in the density structure there (by deepening the poten- 
tial) may have a significant impact on the torque values. We can 
compare our values of the torque with the well known formulae 
for the s pecific torque in a 3 D strictly isothermal disc as pre- 
sented bv lTanaka et 311 ( 120021) 



with 



/r = (1.364 + 0.541 as) 



(18) 



(19) 
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Fig. 9. Radial density distribution in the equator along (p(r, - 
n,0 - 71 12)), i.e. along a ray through the location of the planet 
for all 4 planetary potentials used. 



this horseshoe region will be mixed thoroughly, the torques de- 
cline and settle eventually to their final equilibria, here reached 
after about 40 orbits. This process of phase mixing inside the 
horseshoe region is called torque saturation, and it occurs on 
timescales of the order of the libration time, which is given by 



Tub 



4ap 
3x, 



(17) 



where aj_ denotes the radial gradient of the surface density 
through S oc r""^ . For our standard parameter this formula gives 
about r^"' = -2.5 10"^fl^Q^„ which is, in absolute value, about a 
factor 1.4 - 2.2 times larger than our results. We note however, 
that eq. (fTSl l has been derived for constant temperature, inviscid 
disc. The influ ence of viscosity on the torque has been studied by 
iMassetl (120021) . who found that a reduction of the viscosity from 
our used value of 10"^ to zero will fully saturate the (vortensity- 
related) corotation torque, leading easily to a reduction of the to- 
tal torque by a factor of two. Additional simulations with much 
smaller viscosity (not shown here) indicate indeed that then the 
equilibrium torque is in good agreement with eq. ( fTSl l. 

It seems at first surprising and unpleasing that the torques 
depend so much on the treatment of the planetary potential. 
However, an e-potential has an influence far beyond the Roche- 
radius of the planet and certainly will change the torques act- 
ing on the planet. Here the corotation torques are affected most 
prominently and become more and more positive as the smooth- 
ing l ength is lowered (see also iPaardekooper & Papaloizoul 
|2009|) . Nevertheless, in two-dimensional simulations it has be- 
come customary to rely on e-potentials for the purpose to 
take into account the finite thickness of the disc. In a three- 
dimensional context, the more localised cubic -potential with its 
finite region of influence may be more realistic. But for the 
isothermal case the increased potential depth leads to a very 
large accumulation of mass, as seen in Fig. |9] In such a case 
it will be very difficult to achieve convergence. In the more re- 
alistic radiative case the situation is eased somewhat through a 
temperature increase near the planet, as outlined below. 
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Fig. 10. Specific torque (in units of fl^Qp) acting on the planet 
using 4 different smoothings for the planetary potential in the 
isothermal case with H/r - 0.05. Top: Evolution of total torque 
with time. Bottom: Radial variation of the specific torque den- 
sity r(r) at f = 80 orbits. 



Fig. 11. Specific torque acting on the planet using different grid 
resolutions for the isothermal case with Hjr - 0.037. In all cases 
the cubic potential with rj,„ = 0.5 has been used. Top: Evolution 
of total torque with time. Bottom: Radial variation of the specific 
torque density, at f = 80 orbits for the standard and medium 
resolution, and at 60 orbits for the high resolution. 



To check numerical convergence we performed additional 
runs using a larger number of gridcells. In Fig. [TT]we display 
the total torque versus time and the radial torque density r(r) for 
different grid resolutions, again for the isothermal disc case. In 
contrast to the previous plot we use here a slightly cooler disc 
with Hjr - 0.037, as this matches more closely the results from 
the fully radiative calculations presented below. In this case the 
initial unsaturated torques reach even positive values due to the 
smaller thickness of the disc. The grid resolution seems to be 
sufficient for resolving the structures near the Roche-lobe. For 
the displayed T{r) distribution in the lower panel both cases are 
very similar and the higher resolution case is a bit smoother 

4.3. Fully radiative discs 

The simulations are started from the radiative disc in equilibrium 
as described above, and are continued with an embedded planet 
of 20Mearth- The obtained equilibrium configuration for the sur- 
face density and midplane-temperature is displayed in Fig. [12] 
after an evolutionary time of 100 orbits. As in the isothermal 
case the density within the Roche lobe of the planet is strongly 
enhanced for the deeper potentials, displayed are the two ex- 
treme cases of our different potentials. Comparing with the cor- 
responding density maps of the isothermal case |8] one can also 
observe slightly smaller opening angles of the spiral arms in the 
radiative case. For identical H/r the sound speed would be sfy 
times larger in the radiative case leading to a larger opening an- 



gle. Here, the effect is overcompensated by the reduced temper- 
ature (lower thickness) in the radiative case. A different opening 
angle of the spiral arms will affect the corresponding Lindblad 
torques acting on the planet. 

At the same time, a slight density enhancement is visible 
'ahead' of the planet {ip > 180°) at a slightly smaller radius 
(r<l). This feature that is not visible in the isothermal case is 
caused by including the thermodynamics of the disc. Let us con- 
sider an adiabatic situation just after the planet has been in- 
serted into the disc, and follow material on its horseshoe orbit 
(in the co-rotating frame) as it makes a turn from the outer disc 
(r > 1,1^ > 180°) to the inner (r < 1). The radial temperature 
and density gradient imply for our ideal gas law a gradient in the 
entropy function S in the disc through 




As shown above, in our simulations we find for the surface den- 
sity E oc r"'/- as due to the assumption of a constant viscos- 
ity, and the midplane temperature follows T oc r i-^. Due to 
this gradient in 5 a parcel (coming from outside) has in our 
case a smaller entropy than the inner disc which it is enter- 
ing. Now the entropy remains constant on its path, due to the 
adiabatic assumption. Additionally, dynamical equilibrium re- 
quires that the pressure of the parcel does not change signifi- 
cantly upon its turn, and entropy conservation then implies that 
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the density has to increase. At the same time the density 'be- 
hind' the planet {ip < 180° and r>l) will be lowered by sim- 
ilar reasoning. Both components produce a positive contribu- 
tion to this entropy-related corotation torque that acts on the 
planet, and which adds to the negative Lindblad torque and 
the positive vortensity-related corotation torque. In truly adia- 
batic discs thi s effect will disappear a fter a few libration times 
dBaruteau & Masset 2008; Kiev & Crida 2008.) because the ma- 
terial, being within the horseshoe region, will start interacting 
with itself, and the density and entropy will be smeared out due 
to the mixing, leading to the described torque saturation pro- 
cess. Adding radiative diffusion will prevent this and keep the 
entropy-related torques unsaturated, and a non-zero viscosity is 
also required. 

In this fully radiative case the temperature within the Roche 
radius of the planet has also increased substantially due to com- 
pressional heating of the gas (lower two panels in Fig. [T2b . In 
addition, the temperature in the spiral arms is increased as well 
due to shock heating. 

The density and temperature runs in the disc midplane along 
a radial line at = 180° cutting through the planet are dis- 
played in Fig. [13] As in the isothermal case, deeper potentials 
lead to higher densities within the Roche lobe. The increase is 
somewhat lower because now the temperature is higher as well 
due to the compression of the material. The larger pressure low- 
ers the density in comparison to the isothermal case. Interesting 
is that the maximum temperature is substantially higher than in 
the ambient disc even for this very low mass planet of 20Mearth- 
Considering accretion onto the planet the increase in temperature 
might be even stronger due to the expected accretion luminosity. 

4.4. Torque analysis for the radiative case 

In the upper panel of Fig.[T4]we display the time evolution of the 
total specific torque acting on a 20 Meanh planet for the full radia- 
tive case. In contrast to the isothermal situation, all four poten- 
tials result now in a positive total torque acting on the planet. As 
in the previous isothermal runs, the torques reach their maximum 
shortly after the onset of the simulations (between f a; 10 - 20) 
and then settle toward their final value. In the corresponding 
isothermal case with Hjr - 0.037 the difference between the 
initial positive unsaturated torque and the final saturated value 
has been very pronounced (see Fig. fTTT i. In contrast, in this fully 
radiative case the inclusion of energy diffusion and the subse- 
quent radiative cooling of the disc will prevent saturation of the 
entropy-related corotation torque, resulting in a positive equi- 
librium torque. Very similar results have been found previously 
in the fully radiative regime in 2D simulations (Kiev & Crid^ 
|2008|) . It is important to notice, that the two cubic -potentials 
(which are more realistic in the 3D case) yield very similar re- 
sults. The more unrealistic 6-potentials show rather strong devi- 
ations because, due to their extended smoothing of the potential, 
they tend to weaken in particular the corotation torques which 
originate in the close vicinity of the planet. In the lower panel of 
Fig. [14] the radial torque distribution is displayed for the same 
4 potentials. In comparison to the corresponding plot for the 
isothermal H/r - 0.037 (see Fig. [TTT i we notice that the regu- 
lar Lindblad part is slightly reduced in the radiative case due to 
the larger sound speed. 

Additionally, clearly seen is the additional positive contribu- 
tion just inside r = 1 which appears to be responsible for the 
torque reversal. This feature is caused by an asymmetric distri- 
bution of the density in the very vicinity of the planet, see also 
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Fig. 12. Surface density (upper two panels) and temperature in 
the equatorial plane (lower two) for fully radiative simulations 
at 100 planetary orbits. Displayed are results for the shallowest 
and deepest potential. Upper panels refer to the e-potential with 
- 0.8, and lower to the cubic -potential with r„„ = 0.5, re- 
spectively. 
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Fig. 13. Radial density (top) and temperature (bottom) distribu- 
tion in the equator along a ray through the location of the planet 
for all 4 planetary potentials used, for the fully radiative case. 



Fig. 14. Specific torques acting on a 20 Meaith planet for different 
numerical potentials in the fully radiative case. Top; Evolution 
of total torque with time. Bottom: Radial variation of the specific 
torque for f = 80 orbits. 



Fig. [17] below. It pulls the planet gravitationally ahead, increas- 
ing its angular momentum, leading to a positive torque. Above 
we argued that this eff'ect may be due to the entropy-related coro- 
tation torque of material moving on horseshoe orbits (see also 
iBaruteau & Masseil2008l) . Due to the symmetry of the problem 
one might expect a similar feature caused by the material mov- 
ing from inside out. However, there is no sign of this present in 
the lower panel of Fig. [14] To analyse this asymmetry, we per- 
formed additional simulations varying the grid resolution and 
disc thermodynamics. That the feature is not caused by lack of 
numerical resolution is demonstrated in Fig. [15] where results 
obtained with two different grids are displayed. Both models 
show the same characteristic torque enhancement just inside the 
planet. In Fig. [16] we compare the radial torque density of the 
isothermal and a new adiabatic model for Hjr — 0.037 with the 
fully radiative model, all for the cubic potential with Vsm - 0.5 
at intermediate resolution. For the adiabatic case we show r(r) 
at two different times. The first at r = 10 when the torques are 
unsaturated, and the second at f = 80 after saturation has oc- 
curred. Please note, that the isothermal and adiabatic models 
start from the same initial conditions (locally isothermal), while 
the radiative model starts from the radiative equilibrium without 
the planet. While the adiabatic model at f = 10 shows signs of 
the enhanced torque just inside the planet, there is no sign of 
a similar feature at a radius just outside of the planet. Hence, 
this asymmetry of the entropy-related corotation torque is visi- 
ble in both the adiabatic and radiative case. Outside of the planet 
the adiabatic model and the radiative agree very well as Hj r is 
similar, while the isothermal model deviates due to the different 



sound speed. Whether the location of the maximum in the torque 
density is identical in the radiative and adiabatic case is hard to 
say from these simulations because in 3D adiabatic runs the peak 
appears to be substantially broad er with respect to correspon ding 
2D cases. It has been argued bv IBaruteau & Massetl ( |2008|) that 
it should occur exactly at the corotation radius, which is shifted 
(very slightly) from the planet's location due to the pressure gra- 
dient in the disc. In our radiative simulations it seems that the 
maximum is slightly shifted inwards, an effect which may be 
caused by adding radiative diffusion to the models and consider 
discs in equilibrium. An issue that certainly needs further inves- 
tigation. 

In Fig. [T7] we present additional results of supporting 2D 
simulations for the fully radiative case, as shown for lower reso- 
lution in Kley & Crida (2008). All physical parameter are iden- 
tical to our 3D fully radiative case. The top panel shows the sur- 
face density distribution next to the planet. Seen is the density 
enhancement just inside and ahead of the planet, and some in- 
dication for a lowering outside and behind. Please note, that the 
planet moves counter-clockwise into the positive ^-direction, i.e. 
upward in Fig. [17] In the lower panel we display for each grid- 
cell the net torque (f ^) acting on the planet. It is constructed by 
adding each cell's individual contribution to the torque and that 
of the symmetric cell with respect to the planet location, i.e. 

f ^ = + [r(r/, 4> planet + 4>j) + ^H, 4> planet ' <t> j)] ■ (20) 

Hence, in absolute values the bottom half of the plot (with 
4> < ippianet) rcscmblcs exactly the top half {(f) > (ppianei)- The 
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Fig. 15. Specific torque acting on the planet using different grid 
resolutions for the fully radiative case. In all cases the cubic po- 
tential with r,,,, = 0.5 has been used. Top: Evolution of total 
torque with time. Bottom: Radial variation of the specific torque 
density at f = 50 orbits. 



colours are chosen such that blue refers to negative Vfj and yel- 
low/red to positive values. The signs of Vfj are chosen such that 
the top left and lower right quadrant have the correct sign (+) 
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Fig. 16. Radial variation of the specific torque acting on the 
planet using different thermodynamical disc models. In all cases 
the cubic potential with rs„, = 0.5 and standard resolution have 
been used. The adiabatic model at f - 10 corresponds to the 
time where the corresponding total torque has its maximum. The 
models at J = 80 have all reached their equilibria. 
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Fig. 17. Results of a 2D fully radiative model using a resolution 
of 5 12 X 1536 gridcells. The gray dot indicates the location of the 
planet, and the curve its Roche lobe. The solid black lines show 
the streamlines. Top: Perturbed surface density with respect to 
the case without an embedded planet. The values are scaled as 
S'^^. Bottom: The net torque acting on the planet caused by the 
mass in each individual gridcell using the prescribed smoothed 
torque cutoff function with b = 0.8, see explanation in text. The 
values are scaled as (f^)'^^. 



and the other two are just reversed. Due to this redundancy in 
the plot, only the upper left and lower right quadrant should be 
taken into account to estimate global effects. The mirroring pro- 
cess at the (f) - TT-line (with the mirroring of the colour scale) 
allows an easy evaluation and comparison of the individual con- 
tributions. One can notice that the net torque will be positive due 
to excess material just ahead and inside of the planet. From the 
plot it is also clear that there exists indeed an asymmetry of the 
torques induced by horseshoe material coming from outside-in 
versus material turning inside-out. In the figure, there is only a 
weak indication of a marginal positive contribution just below 
the planet. In additional simulations for purely adiabatic discs 
with different (positive and negative) entropy gradients which 
have either constant density or temperature it has become ap- 
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Fig. 18. Perturbed entropy (Top) and perturbed density 
(Bottom) for the 2D fully radiative equilibrium model using a 
resolution of 512 x 1536 gridcells. The values are scaled as 2'^"^ 
and S respectively. 



parent that the asymmetry is caused by the entropy gradient. In 
the case of a negative entropy gradient (as in our fully radiative 
model) the positive excess torque comes from inside/ahead the 
planet, while for a positive gradient the negative excess torque 
comes from outside/beh ind the planet. Whether t he maximum 
of r(r) lies at corotation dBaruteau & Masse3l2008l) or is slightly 
shifted when radiative effects are considered may deserve further 
studies. 

In Fig. [18] we show the perturbed entropy and density in the 
2D fully radiative model in equilibrium, for a larger domain. 
Caused by the flow in the horseshoe region, there is an entropy 
minimum for larger <p inside of the planet, and a maximum for 
lower (p outside, for r + rp. Both lie inside the horseshoe region 
and close to the separatrix. The overall entropy distrib ution is 
very similar to that found bv lBaruteau & Massetl (|2008|) for adi- 
abatic discs shortly after the insertion of the planet. Due to the 
included radiative diffusion this effect does not saturate in our 
case, and we clearly support their findings even for the long term 
evolution. The disturbed entropy shows in fact a slight asymme- 



try (in amplitude) with respect to the planet, that reflects back on 
to to the density distribution (bottom panel). 

4.5. Planets with different planetary masses 

Following the results obtained in the previous sections we adopt 
now the cubic r,,,, = 0.5 planetary potential assuming that it is 
closest to reality, and study the effects of planets with various 
masses in fully radiative discs. Starting from the 2D radiative 
equilibrium state (see section 14.11 ) we now place planets with 
masses ranging from 5 up to 100 Mganh in the initially axisym- 
metric 3D disc. The numerical parameters for these simulations 
are identical to those discussed above. 

In recent 2D simulations of radiative discs with embedded 
low mass planets the torque acting on the planet depends on the 
planetary mass in such a way that for planets with a size smaller 
than about 40 earth masses the total torq ue is positive implying 
outward migration dKlev & Cridal 120081) . For larger masses the 
forming gap reduces the contribution of the corotation torques, 
and the results of the radiative simulations approach those of the 
fixed temperature (locally isothermal) runs. Our 3D simulations 
show indeed very similar results for planets in this mass regime, 
see Fig. [19] Planets in the isothermal regime migrate inward 
with a torque proportional to the planet mass squared, as pre- 
dicted for low mass planets undergoing Type-l-Migration. Note, 
that we use in these models the temperature distribution for a 
fixed H/r - 0.037. The values for the three lowest mass planets 
(with 5,10, 15Meaiih) are not as accurate due to the insufficient 
grid resolution, remember the 'kink' in Fig. [TT]which refers to 
20Meaith at standard resolution. For the fully radiative disc the 
planets up to about 33Mearth experience a positive torque, while 
larger mass planets migrate inward, due to the negative torque 
acting on them. 

When comparing t he 3D torques to the corresponding 2D 
values as obtained by iKlev & Cridal (l2008l) for the same disc 
mass and opacity law, we note two differences: /) The absolute 
magnitudes of the torques in the radiative case are enhanced in 
the 3D simulations with respect to the corresponding 2D results, 
resulting in even faster outward migration of the planets. This 
result can be explained by the reduced temperature (i.e. verti- 
cal thickness) of the 3D disc with respect to the 2D counterpart 
(cf. Fig. [3]), as a re duction in H typically increases the torques 
( Tanaka et al.ll2002l) . The turnover mass from positive to neg- 
ative torques is reduced in the 3D simulations. This effect is 
caused again by the reduced disc thickness, as now the onset 
of gap formation {Rfj ~ H) occurs for smaller planetary masses. 
The different form of the potential and the softening length may 
also play a role in explaining some of the differences observed 
between the 3D and 2D results. 

Finally, in Fig.|20]we show that the position of the maximum 
of the radial torque density is independent of the planet mass, 
and is therefore a result of the underlying disc physics. 



5. Summary 

We have investigated the migration of planets in discs using 
fully three-dimensional numerical simulations including radia- 
tive transport using the code NIRVANA. For this purpose we have 
presented and described our implementation of implicit radia- 
tive transport in the flux-limited diffusion approximation, and 
secondly our new FARGO-implementation in full 3D. 

Before embedding the planets we studied the evolution of 
axisymmetric, radiative accretion discs in 2D. Starting with an 
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Fig. 19. Specific torques acting on planets of different masses in 
the fully radiative (blue crosses) and isothermal (red plus signs) 
regime. Note, that the isothermal models are run for a fixed 
H/r - 0.037. All torque values are displayed at a time when 
the equihbrium has been reached. 
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Fig. 20. Radial torque distribution in equilibrium for various 
planet masses. The vertical dotted line indicates the location of 
the maximum. 



isothermal disc model having a fixed H/r - 0.05, we find that for 
our physical disc parameter the inclusion of radiative transport 
yields discs that are thinner {H/r = 0.037 at r = 1). We note that 
in the isothermal case the disc thickness is a chosen input param- 
eter, while in the fully radiative situation it depends on the local 
surface density and the chosen viscosity and opacity. Interesting 
is here the direct comparison to the equiva lent 2D m odels using 
the same viscosity, opacity and disc mass ( iKlev & Crida 2008i) . 
as displayed in Fig. [3] Here, our new 3D disc yields smaller tem- 
peratures (by a factor of 0.6-0.7) than the 2D runs. Since the 2D 
simulations have to work with vertically averaged quantities, it 
will be interesting whether it might be possible to adjust those 
as to yield results in better agreement to our 3D results. 

Concerning planetary migration we have confirmed the oc- 
currence of outward migration for planetary cores in radiative 
discs. As noticed in previous research, the effect is driven by a ra- 
dial entropy gradient across the horseshoe region in the disc, that 
is maintained by radiative diffusion. Our results show that plan- 
ets below the turnover mass of about m ^ 33Mearth migrate out- 
ward while larger masses drift inward. The reduced temperature 
in the 3D versus 2D runs has direct influence on the magnitude of 
the resulting torques acting on the planet. As the disc is thinner 
in 3D the resulting torques, corotation as well as Lindblad, are 



also stronger. The turnover mass from outward to inward migra- 
tion is slightly reduced as well for the 3D disc since the smaller 
vertical thickness allows for gap opening at lower planet masses. 
Due to the reduced temperature in the 3D case, the spiral waves 
have a slightly smaller opening angle compared to the isothermal 
case. 

Another interesting, partly numerical, issue that we have ad- 
dressed concerns the influence that the smoothing of the plane- 
tary potential has on the density structure in the vicinity of the 
planet and the speed of migration. We have compared the stan- 
dard e-potentials in contrast with so-called cubic-potentials. The 
results indicate, that a deeper (cubic) potential results in a higher 
density inside the planet's Roche lobe for isothermal, as well as 
radiative discs. Since the potential depth influences the density 
in the immediate vicinity of the planet the resulting torques show 
some dependence on the chosen smoothing length. We note that 
for the more realistic cubic -potential, changing the smoothing 
from 0.8 to 0.5 does not alter the results for the fully radiative 
simulations considerably, and runs at different numerical resolu- 
tions have indicated numerical convergence. Hence, we believe 
that this value is suitable for performing planet disc simulations 
in 3D. The usage of an e-potential cannot be recommended in 
3D. Outside of the planet's Roche lobe one can hardly notice 
a difference in the density structure. Since the cubic -potential 
agrees with the true planetary potential outside of the smoothing- 
length rsm, it is of course desirable to choose this transition ra- 
dius as small as possible, but the achievable numerical resolution 
always will set a lower limit. We found rsm = 0.5 a suitable value 
for our grid resolution and used that in our parameter studies for 
different planet masses. We note, that independent of the chosen 
form of the potential the outward migration of planetary cores 
seems to be a robust result for radiative discs. As the magnitude 
(and direction) of the effect depends on the viscosity and opac- 
ity, further studies, investigating different radial locations in the 
disc, will be very interesting. 

The outward migration of planet embryos with several earth 
masses certainly represents a solution to the too rapid inward mi- 
gration found in this mass regime of classical type-1 migration. 
Growing planets can spend more time in the outer disc regions 
and move then later via type-II migration towards the star On 
the other hand it may be difficult to reconcile this finding with 
the presence of the discovered Neptune-mass planets that reside 
closer to the central star (a x 0.1 AU), but still too far away to be 
ablated by stellar irradiation. 

As seen in our simulations, parts of the disc can display con- 
vection due to the form of the opacity law used. It will be cer- 
tainly interesting in the future to analyse what influence the con- 
vective motions have on the migration properties of the embed- 
ded protoplanets. Additionally, fully radiative 3D-MHD simula- 
tions are definitely required to judge the efficiency of this process 
in turbulent discs. 
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Appendix A: Fargo algorithm 

Multi-dimensional simulations of accretion discs that include 
the y?-direction typically suffer from severe timestep limitations. 
This is due to the fact that the azimuthal velocity is Keplerian 
and falls off with radius. Hence, the innermost rings determine 
the maximum timestep allowed even though the region of in- 
terest lies much further out. One suggestion to resolve this is- 
sue is given by the FARGO algorithm which stands for "Fast 
Advection in Rotating Gaseous Objects" (Masset 2000a). It has 
originally been developed for 2D disc simulations in a cylin- 
drical coordina te system, for details of the implementation see 
iMassetl (l2000allH) . Here, we briefly describe our extension to 
three spatial dimensions in spherical polar coordinates. 

The basic method relies on a directional splitting of the ad- 
vection part, where first the radial and meridional (in direction) 
advection are performe d in the s tandard way. To calculate the az- 
imuthal part we follow lMassetl ( l2000al) and split the angular ve- 
locity into three parts: From the angular velocity of each grid cell 
^i,j,k - iuip)i,j.k/ fi first an average angular velocity is calcu- 
lated for each radial ring ;, which is obtained here by averaging 
over the azimuthal (index k) and vertical (index j) direction 



1 



Z 



(A.l) 



where the summation runs of all azimuthal and meridional grid- 
cells, and N^,Ng denote the number of these gridcells, respec- 
tively. We note, that the summation over the vertical direction 
with index j is not required at this point. In our case, for a thin 
disc where the angular velocity does not vary much with height, 
the vertical averaging simplifies matters somewhat. From this, 
one calculates an integer-valued shift quantity 



Hi = Nint (w,Af/A^) , 



(A.2) 



where Nint denotes the nearest integer function. This corre- 
sponds to a transport by the angular 'shift velocity' 



(A.3) 



Then we calculate the constant residual velocity of each ring 

(A.4) 

and finally the residual velocity for each individual gridcell 

Rewritten, we find for the angular velocity the following expres- 
sion 



(A.6) 



The advection algorithm in the i/j-direction proceeds now in three 
steps. In the first two steps all quantities are advected using the 
standard advection routine with the transport velocities oj':'' and 
wjy^ and then all quantities are shifted by the integer values «, 

in each ring / which corresponds to a transport velocity oj^" . 
Using this splitting, the transport velocities in the advection part 
are given by the two residual velocities w" and (^"J;^, which are 
typically much smaller than ojij^k- Hence, the time step limita- 
tion for the azimuthal direction is determined by the local varia- 
tion from the mean azimuthal flow in the disc which is typically 
much smaller than the Keplerian value. In our case of a 3D disc 
the time step criterion is first given by the normal CFL-criterion 



as presented for example in I Stone & NormanI d 1 9921) where the 
angular velocity w/ ^^ is just replaced by the residual cell val- 
ues wjy^ and wj''. This change provides the major reduction in 
the transport velocity and a substantial increase in the time step 
size. An additional time step limitation is given by the require- 
ment that the shift should not disconnect two ne ighbouring gri d 
cells in the radial and in the meridional direction (lMassetl2000al) . 
Here, this additional limit on the time step reads 



Afs/jfor = 0.5 min 



I'^i.j.k - <JJi-l.i,k\ \<J^i.j,k - <^i.j-l,k\ 



(A.7) 



The second restriction is only necessary in the case, where the 
above vertical averaging in Eq. lA.ll has not been performed. The 
sequencing of the advection sweeps in the FARGO algorithm has 
to be such that the azimuthal sweep comes at the end, hence in 
our simulations we use always: radial, meridional, and finally 
azimuthal. 

In a staggered mesh code such as our NIRVANA code, that 
is essentially based on the ZEUS method, an additional com- 
plication arises in the straightforward application of the FARGO 
method, due to the fact that the velocity variables are located at 
the cell interfaces and not at the centres. Hence, the correspond- 
ing 'momentum cells' for the radial and meridional momenta 
(pui- and prug) are shifted with respect to the standard density 
cells by half a gridcell in the radial or meridional direction, re- 
spectively. To apply the FARGO method one has first to split all 
the momentum cells in two halves, use the algorithm outlined 
above on each of the halves, and then combine them again after- 
wards to calculate from the updated momenta the new velocities 
on the interfaces. This leads of course to an overhead in the simu- 
lation cost which is counterbalanced however by the much larger 
time step. 

Appendix B: Radiative transport 

Here, we outline briefly the method to solve the flux-limited dif- 
fusion equation in 3D. Radiative transport is treated as a sub- 
step of the integration procedure. In equilibrium viscous heating 
is balancing radiative diffusion and to ensure this also numeri- 
cally we incorporate the dissipation into this sub-step. Using the 
appropriate parts of the energy equation ^ and the flux ^ we 
obtain a diffusion equation for the gas temperature. 



dT I ^ 

— ^—[V-DVT + e+] 

Ot CyP 

where the diffusion coefficient is given by 



D 



AcAcirT^ 

p{K + 0-) ' 



(B.l) 



(B.2) 



and denotes the viscous dissipation that is added to the sys- 
tem. The flux-limiter A depends on the local physical state of the 
gas and approaches A - 1/3 in the optically thick parts and re- 
duces the flux to F - ccirT'^ in the optically thin parts. Here we 
use an expressions for A as given in lKlevI d 1989b . 

A straight forward finite difference form of Equation dB.ll ) 
in Cartesian Coordinates is given by 



•,j,k i,j,k 

At 



1 



(CYp)i,j,k 



Ax Ax ■•j-' 



(B.3) 

Ti.j.k - Tj-ijA 
Ax I 
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^ J_ljyy Tjj+ik - Tjj^k Tjj^k - Ti,j-\,k 

l^yy i.j+U ij.k 

1 I ^- Ti,i,k+i ~ Tiik Ti jk-Ti n-\ 

+ £)~ - - £)~ — — 

In orthogonal curvilinear coordinates additional geometry terms 
have to be added in the above equation. Here D^j^. denotes 




i+i,j 



At < min 

i.j.k 



^ ^i.j,k 



where the superscript n + 1 has been omitted on the left hand side. 
The coefficients AJ'^.^ to C^. j^, can be obtained straightforwardly 
from Eq. (IB.3l l. The right hand side is given by 



] 



(B.4) 



and so forth. 

The grid structure from which Eq. ( IB.3l l follows is outlined 
in the two-dimensional case in Fig lB.ll One has to keep in mind 
that the temperature (being a scalar) is defined in the cell centre, 
while the values of ., are defined at the cell interfaces. 



Fig.B.l. The grid structure used in a staggered grid code for 
a two-dimensional example. Shown is the cell i,] where the 
coordinates of the cell centre [jc^yj.] is given by [l/2(jt:,-,y + 

^I'^iyij + yi.i+i)]- The temperature Tij is located at the 
cell centre where also the diffusion coefficient D is defined. The 
averaged diffusion coefficients D are defined at the cell inter- 
faces. 



In Eq. (IB. 3b no time levels are specified on the right hand 
side. For explicit differencing the time level should be f" such 
that the new temperature on the left 7""+' is entirely given by the 
old values T" at time f". This might lead to very small timesteps 
since the timestep limitation is approximately given by 



(B.5) 



where D is given by Dlipc^). 

Hence, often an implicit version of the equation has to be 
used, where all the temperature values Tij^k on the rh.s. are eval- 
uated at the new time or an arithmetic mean between new 
and old times. Even though the diffusion coefficients may de- 
pend on temperature, we always evaluate those at the old time 
f". Otherwise this would lead to a non-linear matrix equation. 

Collection all the terms in eq. (IB.3l l this leads to a linear sys- 
tem of equations with the form 

^/,>,i^''.>.*-i + ^i,j,k'^!di<+i + Bij,kTij,k = Ri,j,k (B.6) 



1" + 



Ql,k 



Written in matrix notation Eq. ( IB.6I ) reads 



(B.7) 



Obviously the matrix M is a sparse matrix with a banded struc- 
ture. Usually M is diagonally dominant but in situations with 
extended optically thin regions this property will be lost. 

Equation ( IB.7I ) can in principle be solved by any linear 
equation package. For simplicity and testing purposes we work 
presently with a standard SOR solver. Using an optimised re- 
laxation parameter w we need about 130 iterations per timestep 
in the initial phase which is far from equilibrium and only 80 
iterations at later times near equilibrium. 

The radiation module of the co de has been tes ted extensively 
in different coordinate systems in iBitschI (l2008h . and we have 
compared our new results on radiative viscous discs obtained 
with the 3D version of NIRVANA in detail with those of the exist- 
ing 2D code RH2D. 
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